Projection method for rapid ab initio calculations of metals 
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I. INTRODUCTION 

The key point in the construction of fast, O (TV) nu- 
merical methods is the gap in the excitation spectrum 
above the ground state which renders the interactions 
short ranged. Systems without gap display infrared sin- 
gular, long range interactions which slow down the con- 
vergence of the numerical algorithms. It seems natural 
to seek different strategies to deal with the short and the 
long range quantum fluctuations. In particular, one may 
use rapid numerical methods for the short range fluc- 
tuations and treat the more difficult long range sector 
with slower, more sophisticated method. Such a mixed 
numerical algorithm is discussed in this paper. 

The strategy of the renormalization groupi is a natural 
candidate for the construction of such an algorithm. In 
fact, the renormalization group is a systematic method 
to successively eliminate certain degrees of freedom or 
fluctuation modes in such a manner that their impact 
on the dynamics is accumulated in the effective theory 
which is constructed for the remaining degrees of free- 
dom. The algorithm proposed in this paper consists of 
two steps. First, a rapid numerical method is applied for 
the elimination of the short range fluctuations. What is 
left is a dynamical problem of the long range fluctuations 
described by an effective Hamiltonian. This problem is 
dealt with in the second step of the algorithm by the 
diagonalization of the effective Hamiltonian. 

The question of central importance for such a mixed 
method is the relation between the total dimension of the 
Hilbert space and the dimensionality of the linear space 
of the effective theory where the exact diagonalization 
is performed. Let us denote by £ an intrinsic energy 
scale of the system and introduce N > and N < as the 
number of modes with energy superior or inferior of £. 
One may call N > and N < ultraviolet and infrared cut- 
off. Our algorithm will be O [N^A but will slow down as 
A^< is increased. Since N < growths with the size of the 
system in the absence of gap and remains finite when the 
gap is present our algorithm might be useful for systems 
with weak gap or truly gapless models of finite size. The 
numerical efficiency compared to other methods will be 
judged by the prefactor of N > in the required computer 
time so long the system size or the gap is kept fixed. We 



believe that this prefactor will be rather small because 
the modes treated in this step are short ranged. 

The traditional renormalization group methodi con- 
sists of the repeated application of a three step procedure. 
The first step is the blocking, the elimination of certain 
variables from the system. This is usually achieved by 
the lowering of the ultraviolet cut-off, the highest energy 
the fluctuations may reach in the system. The second 
step is the construction of the effective theory for the re- 
maining modes. Finally, in the third step which gave the 
name of the procedure, one performs a rescaling of the 
energy or other scales of the effective theory in order to 
restore the original ultraviolet cut-off. This last step is 
not always necessary. 

There have already been a proposal in the literature 
for a partial implementation of this idea, the so called en- 
ergy space renormalization group-'-, realizing the block- 
ing and the rescaling steps. In order to render this scheme 
systematical one can not be content with the naive elim- 
ination of the unwanted modes, the restriction of the 
Hamiltonian into a subspace, but should realize the sec- 
ond step as well, the accumulation of the effects of the 
excluded directions within the subspace retained. For 
this purpose we use a projection method developed in 
nuclear physics^*^. 

In sec. [n] we expand the density matrix formalism, 
which is the foundement of ab initio algorithms. The 
locality principle and its use in linear-scaling methods 
are presented in sec. IIIII An exemple of such algorithm, 
the so called Fermi Operator Expansion is presented in 
section Hvl In the last section we develop a Numerical 
Renormalization Group method in Hilbert Space around 
the Fermi level and propose an improvement inspired by 
projection method. 



II. DENSITY MATRIX FORMALISM 

In the framework of Density Functional Theory 
(DFT)^*2i£, particularly in the Kohn-Sham scheme 8 , 
rapid ab initio calculation methods allowing linear scal- 
ing or near-linear scaling computation time have been 
developed recentlj*i2iii. Most of the rapid ab initio algo- 
rithm is based on the one-particle reduced density matrix 
p which is assumed to be a projector on the subspace 
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spanned by the low-lying occupied states according to 
the auf bau principle : 



= ^2 f 00,^)1^) (ipi 
i 

N/2 



(i) 



where N denotes the number of electrons and (e*, l^i)) is 
an eigenfunction of the Kohn-Sham Hamiltonian H and 
//3 is defined as the Fermi-Dirac distribution function at 
the inverse temperature (3, 



where a 



for the tight-binding limit and 



oi ~ o-iattice ■ Ae gap for the weak-binding limit. For sys- 
tems with no gap the decay of p is algebraic at zero 



p{r, r 1 ) = k, 
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(10) 



Such an algebraic decay reflects the presence of long 
range correlations and prevents linear-scaling in the nu- 
merical calculations. The electron states tend to be more 
localized for disordered systems and the matrix elements 
of the density matrix decay faster with the distanc e) 1 ?' 2 ? . 



(2) 



The form H = ^ ti\ipi)(ipi\ allows us to write the density 



matrix as 



P = foo,e F {H) 

= Q(e F -H) 



(3) 



being the Heaviside function. The average energy and 
the particle number can be written as 



and 



E B s -^PijHj 



i,3 



(4) 



(5) 



where the density matrix is given in a localized orbital 
basis 



(6) 



and Hij = {4>i\H\<j>j) , S i3 = (4>i\4>j)- 

III. PRINCIPLE OF "NEARSIGHTEDNESS" 

This principle statesi2ii2ii4 that the matrix elements 
of the one-electron density matrix are negligible beyond 
the distance ca where a is the lattice spacing, 



> c => pi 3 ~ 0, 



giving 



E 



(7) 



(8) 



i max(0,i—c)<j<min(N 1 i-\-c) 



The decay of the density matrix p in real-space de- 
pends on the material. For systems with gap the decay 
is exponential 12 i 13 i 14 i 15 i 16 



-a.\r— r | 



(9) 



IV. FERMI OPERATOR EXPANSION 

The polynomial expansion of p in Chebychev polyno- 
mials, the so called Fermi Operator Expansion 2 ^ 2 ^ 3 " is 
an important ingredient of rapid algorithms. 

Chebychev polynomials are defined by the recursion 
formula 



7b Or) 
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(11) 



2xT n+1 (x) - T n (x) 



for — 1 < x < 1. It is easy to see that actually T n {x) — 
cos(n arccosa;). We use the functional form J3J in order 
to fit p with a Chebychev polynomials up to order p, 



Pp 



(12) 



where H s is the dimensionless Hamiltonian scaled and 
shifted into the interval [—1,1], 



Hs 



H - E 
AE 



E = — [min(spec(-ff )) + max(spec(H))] (13) 
AE = — [max(spec(Hj) — mm(spec(H))] 

and f3 s = (3AE, p s = (p - E)/AE. The smallest and 
largest eigenvalue of H can be computed by using the 
Lanczos method which scales linearly with the size of the 
matrix. The projection coefficients 



Q"n {Ps ) ps 



2-8, 



nO 



cos( n9) — - r 



(16 



will be computed numerically by means of Fast Fourier 
Transform (FFT). The particle number conservation fixes 
the value of the Fermi energy level tF found by solving 
Eq. ©. 

The accuracy of Fermi Operator Expansion can be es- 
timates by recalling that one truncates the Chebychev 
polynoms T n of Eq. (|llfl in such a manner that only the 
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matrix elements (T n (H))ij with \i — j\ < c are retained. 
The computation time will be of order pc 2 N = O(N). 
It can be shownSiS^ that the order of the Chebychev 
expansion should be 

p~|(d-l)A = |(d-l)/3A£! (15) 

for the accuracy 10~ d of the expansion coefficients {aj}. 
If the system has an HOMO-LUMO gap AE gap and 
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then 



P 



> 



4{d-l)AElog w d 



3AE, 



(16) 



(17) 
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Since the range of correlations in the density matrix is 
bounded, 



range(p) <p c~ -c(d-l)/3AE (18) 

o 

the correlations grows with the inverse temperature for 
gapless systems and the linear-scaling methods are ren- 
dered inapplicable. 



V. ENERGY RENORMALIZATION GROUP 



density matrices p n4L — fp nifJ ,(H). The zero tempera- 
ture expectation value of an observable A is written as a 
telescopic series 



(A) = Tr( Poo ^A) =^Tr(A n ,^A). 



(19) 



where 



Pn 



Pn-1,M> ^0,i 



P0„ 



Each term in this equation corresponds to a more re- 
stricted energy subspaces centered at the Fermi energy 
level as n is increased. The localization in the energy 
leads to delocalized states in real space in the absence 
of disorder. The ground state is approached by the tele- 
scopic summation by zooming onto the Fermi level and 
the corresponding density matrix projects on more and 
more extended states. 

The order of the Chebychev expansion p is chosen to 
be independent of n and the coefficients obtained by FFT 



arc 



where 



A n ,^ — fp n ,n(H n ) - fp n _ UIJ ,{H n ). 



(20) 



(21) 



We present now a renormalization group method in 
the energy space in order to treat systems with small 
gap. In the original version of this method^ one starts 
with a series of inverse temperatures j3 n — > oo and the 
corresponding density matrices p n which tend to be con- 
centrated around the Fermi level as n — > oo. This alone 
would not represent any improvement as far as the nu- 
merical difficulties of obtaining the density matrices are 
concerned. But the density matrices are constructed in 
decreasing subspaces 7i n D H n +i where Ti. n +i is span by 
the eigenvectors of p n with large eigenvalues. 

This algorithm is modified in order to implement the 
blocking in energy space. First a common chemical po- 
tential is introduced for each temperature which is ad- 
justed at the end of the computation to dial the desired 
particle number. This modification is needed to clear 
the way for the blocking. The Hamiltonians were simply 
truncated in the original algorithm as their subspaces 
were restricted. In order to retain the dynamics of the 
excluded dimensions we employ a method developed in 
Nuclear Physics- which yields an exact, O (N 2 ) algo- 
rithm. 



A. Blocking in the Hilbert pace 

A geometric series of inverse temperatures j3 n = q n (3o 
is introduced for q > 1 together with the corresponding 



ERG 




Fermi level 
Energy 



FIG. 1: Spectral representation of p n and A. 



B. Fixed-point 

The convergence of the telescopic series can be ex- 
pressed as the existence of a fixed point of the blocking 
in the energy space for energy dependent operators. In 
fact, let us suppose that a continuous operator A is com- 
muting with the Hamiltonian H and can be diagonalized 
in a basis of eigenvectors of H . We can then express its 
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expectation value by means of p as 

Tr(A n+1 , M A) = J de A(e)A n+1 ^(e) 

Pn f 7 A f fin 



fln+l 
Jn_ 

(3: 



de Ale 



'A 



'n+l 



A„, M (e)(22) 



n+l 



Tr [ A n ^A 



0n 



n+l 



This expression allows us to rescale the operator A 
around the Fermi- level by the factor f3 n+ i/(3 n and to keep 
A„ iM unchanged. Since A is a continuous operator the it- 
eration of this step obviously leads to a fixed-point, 



Tr(A n+Ul A) - Tr(A n<fl A) -► 



(23) 




FIG. 2: Spectral representation of d and 



as n — > oo. 



C. Projection 

The identification of the subspaces proceeds by the 
construction of the projectors P„ : TL n — > TL n +i- We in- 
troduce first the following pseudo-projectors constructed 
by means of the Chebychev expansion 



C>Pr, 



dfi 



= PnPn,^ ~ Pn,a) 



(24) 



Once the series {G n } is found another set of matrices 
{C n } is formed. The columns of C n are basis vectors 
of Tin+i by means of a heuristic version of the singular 
value decomposition with column pivot ing2iiS&. Since 
the dynamics of the excluded dimensions is retained in 
our case this choice is a less sensitive step of the algorithm 
then in the original version and influences the sparsity of 
the resulting density matrices only. As the next step, the 
overlap matrices S n = C*C n are constructed. Finally, 
the projectors are given as = C i S~ 1 C*. S± can 



H n+ \ in the restricted space with the same dynamics 
around the Fermi level as those of H n , 

1 



H n +i — S n 2 C* I H n + H n Q 



P Qn,HnQ'< 
_ l 
X C n S n 

H n +i Sn 5 C n H n Q n 



1 



-Q n H r 

i / 

(27) 

_ i 



p, — H n Q n H n 

where Q n = 1 — P n . The exclusion of directions from 
the Hilbert space renders the finding of the projection 
of the eigenvectors of the original Hamiltonian into the 
restricted space a nonlinear problem. This complica- 
tion appears as a nonlinear dependence of the eigenvec- 
tor equation in the restricted space on the eigenvalue. 
The energy eigenvalue was replaced by the Fermi level, 
p. in the 'self-energy', the second term on right hand 
side of Eq. I|27|) . The inverse in the right hand side can 
be obtained by the well-known Schultz's or Hotelling's 
method 28 i 29 i 30 as (p — Hn)^ 1 = lim^oo Xj where 



x j — Xj 



.i[21 - (ji - HJXj-i] 



(28) 



actually be obtained as 5, 
of the algorithm^ 



1 



-1/2 



limfc^oo Ak by the help with the initial-guess^ X = (p 



A k = -{ZA k - A k B k A k ) 
B k = ^(W k -B k A k B k ) 



(25) 



-y/a ■ 1, Bq 



-yfa ■ Si and a 



with A 

1 / max | (Si)j k | . The projected Hamiltonian is of the form 



TjERG C 2 TT (~i Q 

•"n+l — ° n C/„-fJ rl < -'n»Jn 



(26) 



ff«)7£(M 

The calculation ends when the dimension of the sub- 
space is sufficiently small for explicit diagonalization. 

One can introduce approximations which render the 
method O (N) . One possibility is the note that Xj of 
Eq. (|28|l converges quadratically and the order of 30 
iterations, a value independent of the system size was al- 
ways sufficient in our numerical test. Another possibility 
is based on the adjustment of the chemical potential at 
the end of the computation. This circumstance allows us 
to make the replacement 



1 



Up to now we have excluded certain directions of the 
Hilbert space which are supposed to be less important 
from the point of view of the ground state dynamics. 
In order to perform the analogue of the Kadanoff- Wilson 
blocking we have to construct an effective Hamiltonian 4,5 



I 

/' 



n ii n •• (29) 

in Eq. 127fl where p will include the 'average' of Q n H n Q n 
within TL n . Such a simplification is more acceptable for 
large n where diuiH n is not too large and the evolution 
is slow. 
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TABLE I: Relative errors of energy computations for different size of systems with /3o — 5, q = 10 and with a Chebychev 
expansion p — 10 



N 


/v IT 1 ™ ^. + 

Hi xact 






EttOT tp 




Et'TOTp — 


256 


93.39 


95.44 


91.58 


2.20 


1.94 


0.26 


384 


139.90 


142.95 


139.09 


2.18 


0.58 


1.60 


512 


186.41 


190.47 


186.61 


2.18 


0.10 


2.07 


640 


232 93 


237.98 


234.12 


2.17 


0.51 


1.66 


768 


279.44 


285 50 


281.63 


2.17 


0.79 


1.38 


896 


325 95 


333.01 


329.15 


2.17 


0.98 


1.19 


1 094 


^79 zL7 


OOU. 


o / u.uu 


9 1 fi 


1 15 
J — LO 




1152 


418.98 


428.04 


424.17 


2.16 


1.24 


0.92 


1280 


465.49 


475.55 


471.69 


2.16 


1.33 


0.83 


1408 


512.00 


523.06 


519.20 


2.16 


1.41 


0.75 


1536 


558.52 


570.58 


566.72 


2.16 


1.47 


0.69 


1664 


605.03 


618.09 


614.23 


2.16 


1.52 


0.64 


1792 


651.54 


665.60 


661.74 


2.16 


1.57 


0.59 


1920 


698.05 


713.12 


709.26 


2.16 


1.60 


0.55 


2048 


744.57 


760.63 


756.77 


2.16 


1.64 


0.52 



D. Numerical test 

We considered a lattice of 2N sites in one dimen- 
sion with nearest neighbor interaction described by the 
HamiltonianSi 

H = 2j2af ai - a t*j (30) 

i <i,j> 

at half filling. Being the simplest model for the conduct- 
ing band electrons the matrix elements of the density 
matrix, computed in the appendix for half-filling, show 
metal-like decrease with the distance. 
Table U shows the results of energy calculations with the 
algorithm of Eq. i|27|) for different sizes. It has been 
reported 2 *^ that the CPU time of the ERG method scales 
as N In 2 N. The computation of Eq. Q27[) which was done 
by applying the approximation lf^|l does not change this 
result since it contains matrix multiplications only. 



VI. CONCLUSION 

A new application of the renormalization group 
method is presented in this work. This method is de- 
signed to retain the dynamics of modes excluded from 
the computation and was developed for the path integral. 
But it is an ideal tool to improve systematically the trun- 
cations of the Hilbert space committed in the operator 
formalism, too. As an example the improvement of the 
Energy Renormalization Group was presented. Here the 
Kadanoff- Wilson blocking is performed in energy space 
and the effects of the directions of the Hilbert space lost 
by the truncation is retained. Therefore the dimension of 
the linear space is reduced but the physics which can be 
described by states within the reduced space remained 
the same. As long as the ground state and the low ly- 
ing excitations are kept in the linear spaces constructed 



in this sequence the salient features of the model can be 
described in a systematical and more economical manner. 

The elimination of dimensions makes the eigenvalue 
equation nonlinear in the eigenvalues, an effect which 
is well known in many-body theory. In fact, say the 
self energy of a particle receives a complicated, energy 
dependent contribution from 'virtual', particle-number 
changing processes which leave from and return to the 
one-particle sector of the Fock space. We employed a 
widely used approximation which becomes exact for the 
ground state and the low lying excitations, the replace- 
ment of the energy eigenvalue by the Fermi level in the 
self-energy. The computational need of the resulting 
method is O (N>) with a prefactor which growth with 
the volume. Nevertheless we find this result remarkable 
since systems with small gap can safely be treated by 
exact diagonalization in a low dimensional subspace. 

We employed a further simplification of the effective 
Hamiltonian in our numerical test. We replaced the part 
of the Hamiltonian which belongs to the eliminated direc- 
tions and appears in the self-energy by a 'mean-operator' 
which is proportional to the identity. This approximation 
is supposed to become exact for the ground state and the 
low lying excitations of a Fermi-liquid. The method is 
O {N> In 2 iV>) when this simplification is used. 

Our method was tested numerically in the case of the 
one dimensional tight binding model. The ground state 
energy improved and a reduction of its error by 25% was 
found compared to the original algorithm for N = 2048. 

The main question left open by the present work is the 
dependence of the computational requirement on N < , the 
physical size of the system, and the explorations of al- 
ternative approximations which ultimately speed up the 
algorithm in this respect. 
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APPENDIX A: ONE-PARTICLE DENSITY 
MATRIX 



If /i and v have same parities, i.e. [i — v ia even then 
jO M y = 0. For v = \i + 2k + 1 



This Appendix contains some details of the computa- 
tion of the density matrix for the tight binding model of 
Eq. $® at half filling. 

The matrix of eigenvectors corresponding to the N low- 
est eigenvalues of H is given by 



Can 



1 



N 



2N+1 



(Al) 



for 1 < /j, < 2N, and 1 < v < N. The reduced density 
matrix 



1 



AN + 2 



-i)» 



7T 2fc + l 



2 2JV+1 



sm 5 



2 2AT+1 



(A4) 



In order to find rate of decrease of /9 M „ we consider the 
limit N — > oo but keep v — \x = 2k fixed, 



Pn,n+2k 



(~l) fe 

2fc7T 



(A5) 



P = CC* 



(A2) 
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